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ABSTRACT 

Compact finite difference schemes for hyperbolic and convection- 
diffusion equations are presented and their relationships to box schemes 
are described. A simple modification of the mesh ratio At/Ax is shown 
to make a previously described non-dissipative scheme for the hyperbolic 
problem dissipative. The dissipative scheme for the convective-diffusion 
equation is formally second order accurate for all values of the local 
cell Reynolds number. Applications to nonlinear problems are described. 
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Introduction 


This paper describes implicit finite difference schemes for two 
closely related classes of mixed initial-boundary value problems in one 
space dimension. 

Part I treats the hyperbolic problem: 


(la) 

U + AU = 

o, 

0 < x < 1, 0 < t < T 

with the 

b x 

initial and boundary conditions 

U = 

n°. 

t = 0 

(lb) 

B 0 U = 

g o* 

O 

II 


b 2 u = 

8 1* 

X = 1. 


Here A is a nonsingular (r x r) matrix with k positive and l 
negative real eigenvalues while Bq is (k*r) of rank k and 

is (&xx) of rank £, k + £ = r. The boundary conditions are 
assumed to be dissipative, i.e. for every vector satisfying the homo- 
geneous boundary conditions then (-1) U’AU > 0 for x = 0,1. 

Note that the differential equation in (1) is the nonconservative 
form of the system U t + F x (U) - 0 where A * gradF. 

Part II treats the scalar convective-diffusion equation (v > 0) 


(2) u, + a u - V u =0, 0 < x < 1, 0<t<T, 

\ / i; x xx 7 * 

which we write as the system 


(2a) 


u . + a u - Vv — 0 
t: x x 


u - v = 0; 
x 


with W = (ujv) 1 , initial and boundary conditions are 
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u = t = 0 

(2b) B Q W = g(J , x = 0 

B 1 W " 8 l» x “ U 


where Bq,B^ are both (1x2). We also assume that homogeneous boundary 

x+1 

conditions result in the inequalities (-1) (au-Vv)u >_ 0, x = 0,1. 

We are interested in demonstrating the following: 

a) the dissipative implicit finite difference schemes described 
below allow the related nonlinear conservation forms of (1) and (2) to be 
treated in their nonconservation forms; numerical evidence for this asser- 
tion is provided below. 

b) for the linear hyperbolic equation dissipation can be restricted 
to affect amplitude modulation while phase errors will be affected only 
by the CFL number; it is thus possible to control pre or post oscillations 
at a discontinuity by the choice of the mesh ratio X = At/Ax. 

c) for the convective-dif fusion equation the scheme is second order 
accurate for all values of the "local cell Reynolds number" ( | a j Ax)/ (2v) . 

A common feature of- the schemes developed for both classes of these 
problems is that they are equivalent to box schemes (Keller [5]) and may 
be solved by the algebraic methods described by Keller in [6] (c.f. 

[2], [4], [8]). Such compact schemes for mixed initial boundary value 
problems avoid boundary extrapolation techniques which are generally 
required to make noncompact schemes algebraically determined and which 
can be an important source of errors. 


The following notation will be employed: 

6 x u i = (u i+i" u i-4 )/Ax ’ 6 t u i = (u i " u i )/At ’ 

n , n , n wo .. n , n +2 , n-J\ . 

P x u i = (u i+J +u i-£ >/2 ’ Vi = (u i +u i )/2 * 



( 3 ) 
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Clearly, 


so that 


5(uv) = y(u)<5(v) + y(v)<5(u) 


J6(u 2 ) = y(u)S(u). 


I. The Hyperbolic Problem 
1.1 A Non-Pis sip at ive Scheme . 

In order to make this paper self-contained as well as to motivate the 
extension to be described below this section reviews a treatment of the 
system U + AU^ = 0 when A is symmetric and constant which has been 
given elsewhere (Rose [11]; also see Wendroff [12], [13]). 

Consider the scheme 


( 1 . 1 ) 


(a) 

(b) 


6 U n + A6 U n 

t 1 XI 


T n , xi 

u t u i ■ Vi' 


- o 


Under the conditions just indicated this scheme provides a convergent approxi- 
mation to (1) as the following sketch of an ’'energy 11 argument shows: multiply 
(I. la) by y and employ (I. lb) to obtain 

M(U?)’U”) + A6 ((u")'U n ) = 0. 
fc ±' x x i l 



= I (U?)'U?Ax, 

L X 1 


Summing on i, setting 
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and noting that the boundary conditions (lb) were assumed to be dissipa- 
tive, there results 

(1.2) || U n || < || U°|| . 

This estimate implies that the solution of (1.1) converges to the solution 
of (1) for all values of the mesh ratio X = At/Ax, Ax -*■ 0. A closer 
examination also shows that the scheme (1.1) is non-dissipative. 

The solution of (1.1) may be obtained by one of the following methods: 

(i) Two-Step Method 

(a) Eliminate uT*’ 5 from (1.1) to obtain 

(I - 3 ) V A)u i+i + R_(A)u ”_ 4 = u"' 2 , 

where 

R ± (*) = l(I±XA). 

The difference equations (1.3) present a two-point algebraic boundary value 
problem which is solvable (Keller [6]) under the boundary conditions described 
in (1) with the initial conditions given by U^~ 2 . 

(b) with U # so determined, explicitly solve (I.la) or (I.lb) for 



(ii) One-Step Method (Box Scheme A) 

n_A 

Instead of employing the explicit step (b) eliminate U. 2 
from (1.1) to obtain 

+R + »)d“_j 


(1.4) 
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Next replace n by n+1 in (1.4) and then eliminate from (1.3) 

and (1.4) to obtain 

(1.5) LUJ 4 * 5 R + (X)U“^ + R_(X)U^ - R_(X)U^ - R + (X)u"_^. = 0. 

This algebraic system is similar to (1.3) and is solvable for by 

the same method. Its relationship to a box scheme is evident from the 
fact that the values occurring in (1.5) are associated with the vertex 
points of a computational cell centered at (i,n+£). The implementation 
of this scheme involves first employing step (a) in (i) . 

(iii) One-Step Method (Box Scheme B) 

Equation (I. lb) is identically satisfied by introducing new 
variables V such that 


20 i4i ' V £f + V i4 


and 


2u n+ i _ yTl+ 1 + 

i i i 


Substitution in (I . 1) then yields a box scheme for values centered at (i,n). 

This form of the box scheme method is essentially that developed by Keller 
[5] and will not be discussed in further detail in this paper. 

Considered separately from (i) the box scheme formulations tend to ob- 
scure the existence of the simple energy estimate (1.2). However, the box 
scheme formulation has advantages for the analysis of other questions. Thus, 

an expansion of the operator LU^"^ immediately results in a truncation error 

2 

estimate which is proporational to (Ax) while an analysis of the amplifica- 
tion matrix in (1.5) also shows the scheme to be non-dissipative and, hence, 
will be of limited use in treating discontinous solutions. The following 
section describes a simple modification of (1. 1) .'which" yields a dissipative system 
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1.2 A Dissipative Scheme 

Instead of (1.1) consider the scheme 


( 1 . 6 ) 


(a) 6 t U^ + A<$ x U" = 0 

(b) y t u" = y x uj + aA6 x u" , 


with a >_ 0, 

Assuming again that A is symmetric and constant, and suppressing 
indices, the energy expression resulting from (1.6) is now 

o = hi {6 t (U'U)+A5 x (U , U)} + a l (6 x AU)'(6 x AU), 


so that 


(1.7) 




Clearly this system is dissipative, i.e., when a > 0 the inequality in 

(1.7) holds unless U n is constant. 

Let 

. . . . . £ - 2a /Ax 

and 

(1.8) X ± *X±e. 

The following solution methods now result: 

(i) Two-Step Method 

a) R + (A + )u" + £ + R_(X + )U“_£ - 

(1.9) 2 

b) 6^ + A5 x U" = 0. 
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(ii) 
(I. 10) 


One-Step Method 

L e u f* 5 V x+ > D Sj + M x+ n£j - R - (X ' )u ?+t - R -< x 'K-i - »• 


Thus (1.9) and (I. 10) result from (1.3) and (1.4) by appropriately 
increasing or decreasing the parameter X. If W? = - U, where U is 

the solution of the differential equation in (1), the estimate of trunca- 
tion error using (I. 10) is 

(At) -1 L £ W^ = 20A 2 U xx + 0(Ax 2 ), 

so that the parameter a = eAx/2 gives rise to an artificial viscosity term 
proportional to U 

XX 


I. 3 Amplification and Phase Error 

Consider the scalar equation u f +au = 0 where a > 0. Initial data 

x 

given by u^ = exp(i0x) are transformed by (1.6) into v = p exp[i(0x-^)] 
while the differential equation carries u^ into v f = exp [i0 (x - aAt) ] . 

Eq. (I. 10) shows that 


(I. 11) 


while 


2 _ 1 + (aX~) 2 

M +2 

1 + (aAV 


< 1 , 


(1.12) = arc tan(aX tan(0/2)). 

Thus | p | < 1 for e > 0, i.e. (1.6) is dissipative for a > 0. 

For aX = 1, = 0. Since tanacj) > atancj), 0 < a < 1, then ij; < 0 

for 0 < aX< 1; similarly iJj > 0 for aX > 1. Thus the wave velocity of 



the wave solution of the differential equation according as the CFL 
number is the difference equation has the same relationship to that 
of less than, equal to, or greater than 1. 

In the non-scalar case an analysis of the amplification matrix 
associated with (I. 10) yields an inequality similar to (I. 11). This 
observation also implies the convergence of the scheme (1.6) since it 
is, clearly, consistent with the differential equation (1) . 


1.4 Numerical Experiments 

The preceding discussion concerned the case A = constant. For 
nonlinear problems, A = A(U) and it is natural to apply (1.6) in 
which the coefficient A is determined by U^. Because (1.6) is 
equivalent to an artificial viscosity method it may be expected that 
the dissipative scheme (1.6) will converge to the physical weak solution 
of the nonlinear conservation system U^. 4- ^ X (U) = 0. 

In one dimension the nonconservation form of the Euler equations 
for inviscid fluid flow is described by 

11 + AU =0, 
t x 

where U = (p,u,p) 1 (p = density, u - velocity., p - pressure) and 


r 

(1.13) A mi 0 

\0 

(Y = 1.4). 

*' 4 The results of several numerical experiments with Riemann problems 
employing the two-step method (1.9) will be presented here. 

Figure 1 illustrates the numerical density profile of a shock 
travelling to the 'right with speed 0.979 which results from the initial 



conditions 
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x < 0 

p = 0.313 
u = 0.3 

p = 0.166 


x > 0 

0.219 

0.0 

0.1 


The Indicated values of p and u on the left and p on the right were 
used to supply boundary conditions. The dissipation factor in (1.9) was 
e = 0^15 and Ax = .01. The value X = 1.04 in Figure la approximates 
the situation in which the average value of the CFL numbers before and 
after the shock was 1. The smoothness of the transition across the shock 
and the fairly accurate tracking of the correct shock position (indicated 
by the vertical line) are evident. Figure lb illustrates the result of 
increasing the CFL number on both sides of the shock (X =1.3) while in 
Figure lc the average CFL number was reduced (X * 0.7). The post shock 
oscillation when X = 1.3 and the preshock oscillation when X = 0.7 
which was evident in the figures may be interpreted as the influence of 
the CFL number on the wave velocity as discussed in 1.3. 

Figure 2 illustrates the p, u, p profiles resulting from the initial 
conditions 


x < 0 
p = 1.0 
u = 0.0 

p = 1.0 


x > 0 
0.125 
0.0 
0.125 


The analytical solution, represented in the figures by the continuous line, 
yields a shock with speed 1.822 and a contact with speed 0.878. The calculated 
values at time = 2.4 with Ax = 0.1, X = 0.6, and a dissipation factor of 
£ = .125 yielded results in fairly good agreement with the exact solution. 

The relative error of the shock speed was estimated to be 3% and the relative 
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error of the total energy was calculated to be 0.3%. The value 
X = 0.6 chosen results in an average CFL number = 1 through the shock 
zone. In this experiment the matrix A in (I. I) was estimated by 
aJ - ilACU^)] + 4[A(uJ_^)]. 

In both of the above experiments the dissipative factor e was 
kept constant for all points. 

II. Convective-Diffusion Equation 
II. 1 Difference Scheme 

A motivation for this scheme will be given in the Appendix* 

Consider the scalar convective-diffusion equation in the form 
described by (2a), i.e. 


(II. 1) 


Let 


u t + - OT x " 0 

u x - y - 0. 

6 = aAx/2v, 


so that | 6 | is the local cell Reynolds number. 
Also let 


(II. 2) 


p (0) = 0 2(0 coth 0 - 1) , 

q(0) = 0p(0).. 


The following approximations to p and a .will be convenient: 

(i) 1 0 1 < 3 : p n, 1/3, q % 0/3 , 

(ii) 1 0 1 >_ 3 : p 1 0 1 \ q 'u sgn 0 - 1/0 . 
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The proposed scheme for treating (II. 1) is: 



a) 

* 11 
O^U . 

+ a6 Uj - v6 v” 

= 0 


t 1 

x i x i 


(II. 3) 

b) 

n 

y u, 
x i 

,Axv 2 r n 
- ( T ) P Vi - 

n 

y t u i 


c) 

* n 
0 u. 

X X 

+ ( T ) q Vi ’ 

y v?. 

X x 


For convergence Ax 0 in which case the coefficients of p and 
q in (II. 3) may be neglected. For | 0 | large, however, the term in- 
volving p and q can provide important corrections to the scheme. 

The difference equations (II. 3) are clearly consistent with the 
differential equation (2). The convergence of the scheme for Ax 0 
when a = constant, is most easily shown by employing the following "energy” 
argument: neglecting the terms in (II. 3) which involve p and q and 
suppressing the indices (i,n), multiply (II. 3a) by li u and then employ 
(II.3b,c) to obtain 

0 = 26 u 2 + Ja 6 u 2 - vS (uv) + v(6, u) 2 . 

t XX X 

Assuming the boundary conditions satisfy the dissipative inequalities 
related to ( 2 ) as described in the introduction, summation over the 
spatial index then yields 

<n.*> ||u n || < ||u°||, 

where 

II u"|| - l (u “) 2 • ta. 

i 

(The inequality in (II. 4) is strict unless u^ is constant). This implies 
the convergence of solutions of (II. 3) to the solution of (2) as Ax -> 0 for 
all values of the ratio X ~ At/Ax. V 



- 12 - 


II. 2 Solution Methods 

The solution method described earlier for hyperbolic systems has 

a direct extention to the equations (II. 3). It will be convenient to 

T 

introduce the vector W = (u,v) and the matrices 


(l±aX) 



±1 


+ (v\ + ~-p) 

Ax,. _ . 
- ~2~(\ + q) 


and 


(1+ aX) 



±<v*-fp> 

- + q) 


The following solution methods result; 


(i) 


Two-Step Method 

By eliminating the term ip (II.3a,b) there results 


(II. 6) 


R X-n + 


R WJ , 




This is an algebraic two-point boundary value prpblem which is solvable 
under the boundary conditions (2b) and the initial values u^ 2 , With 
determined, u^*^ may be obtained explicitly from (II. 3) . 

(ii) One^-Step Method (Box-Scheme) 

By eliminating the term u^ ^ in equations (II. 3a) and (II. 3b) 
and then relabeling the index n the following bpx scheme results: 
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(II.7) Luf* ; - < s + w ?+i +s - H U> ■ °- 

2 

The truncation error, as estimated from this expression, is 0(Ax ) 
for all values of the local cell Reynolds number 1 0 1 . 

II. 3 Numerical Results 

The solution u = cos (x - t)exp (-Vt) of u + u m Vu =0 was 

t X XX 

-2 

computed at t = 4 tt using II. 3 for V = 10 and the values 

Ax = tt/20, tt/40, 7t/80. The norm of the numerical errors as a function 

of X = At /Ax are given in the following table: 


Ax 

X = 

0 , 

.5 

A = 1.0 

A = 2.0 

ir/20 

1. 7 

X 

io“ 3 

0.5 x io -3 

9.3 x io -3 

tt/40 

0-3 

X 

10- 3 

0.27 x io -3 

2.5 x io -3 

tt/80 

0.8 

X 

io -4 

0.07 x 10 -3 

0.6 x 10 -3 


The results confirm the assertion made earlier that II. 3 is second order 
accurate. Similar results were obtained using the norm. 

Figure 3 indicates steady-state boundary layer profiles (t=10) for 
u + au - Vu “0 with boundary conditions u(0) = 1, u(l) = 0 and 

t X XX 

initial condition u(x,0) » (1-x); (II. 3) was employed with Ax *= 1/20 and 
X =* 1. The exact solution is indicated by the solid curves. The numerical 
results indicate fairly close agreement with the exact solution in the 
neighborhood of x = 1. 

2 

Figure 4 describes results for Burger’s equation u + (u /2)x = Vu 

t xx 

with u(x,0) = 1 for x < 0.5 and u(x,0) - 0 for x > 0.5 for 
V = 10 2 , 10 3 at time t = 1.0 using (II. 3) with Ax =* 1/50 and 
X = 1.0; the vertical line indicates the position of the shock for the 


limiting value V « 0. 
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Figure 5 illustrates the solution of Burger's equation at time 
t = 2.0 (steady-state) with boundary condition u(0) - 1, u(l) = -1 
and initial condition u(x,0) = l-2x. The maximum value of the local 
cell Reynolds number was R c = 


112.23. 
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Appendix A 

The Underlying Approximation Method 

The schemes described in this paper have their origin in a common 
approximation method. 

Divide the fundamental solution domain V 0<x<l, 0<t<T 
uniformly into M*N rectangular cells each of area AxAt. If tt” is a 
typical cell with center point (x^,t n ) let a ”+j> T^ -2 denote its 

vertical and horizontal sides as indicated in Figure 6. 



n 


i+i 


Figure 6 


There are thus a total of (M+1)N + (N+1)M sides of which 2(M+N) 
lie on the boundary of V and 2MN - (M + N) lie interior to V. 

n±4 

Consider first the equation U t + AU x =0. If are the average 

values of U on the sides T^ -a , of Tr” then the MN conditions 


i±i 


(A.l) 


6 t u“ + A6 x u" = 0, 


will imply that Gauss* theorem holds on the union of any contiguous set 
of cells • Suppose now that the global solution is approximated by. functions 
which are solutions of the linear differential equation in each cell each of 
which depends, say, upon a parameters, i.e., if are 

linearly independent solutions in a cell, we let 
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(A. 2) 


a 


U - I -'jLp 

i=l X 


i* 


If a = 2 the mixed initial-boundary value problem (1) may be 
approximated as follows (c.f. [10]): set ^ = = ( x ^ -t A); 

then the parameters c^, c ^2 will be determined by any two of the 
four average values U^ -2 and associated with the sides of 

it®. Elimination of the parameters yields two relationships between 
these values one of which is expressed by (A.l) and the other by 


(A. 3) 


^ ’ Vr 


There thus result 2MN conditions for the 2MN + (M+N) average 
values. By imposing the boundary and initial conditions in (1) a 
determined system of equations results. As we have shown through the 
use of an energy argument, when A is symmetric and constant this 
approximation method converges in an norm for smooth solutions. The 

numerical results presented earlier indicate that the dissipative 
scheme based upon a modification of this method provides accurate 
approximations to nonlinear problems as well. 

For the scalar equation u^ + au - vu = 0 similar arguments 

t X XX 

show that, if v * u , then 

x 


(A. 4) 


5 u! + a5 u” - vS v \ 
t i x i x i 


0 , 


is a necessary condition that Gauss' theorem hold in terms of the boundary 
data of cells. Employing (A. 3) with a =3 and the elementary solutions 



±2 = 

i 3 = 


(x- at) 

ax/v 

e 


(A. 5) 


7 



an elimination of parameters yields the scheme (II. 3) which is a determined 
algebraic system under the initial-boundary values given by (2b), Again, 
the energy argument given earlier establishes the convergence of this approxi 
nation scheme when a is constant. 

Both schemes employ an approximation basis which consists of wave 
solutions of the form 

4 >(B,y) = exp[ 3 (x- yt)] , 

where y is the wave velocity. In the hyperbolic case y = A and the 
polynomial basis (l,xl - tA) results by setting 4^ = <f>(0,A) , <j > 2 = 4>g(0,A). 

For the convective-diffusion equation the dispersion relationship 
y = a - v3 holds. In addition to the polynomial solutions (l,x-at) 
the function exp(ax/v) provides another linearly independent solution 
for which y = 0 when 3 = a/v. The basis (1, exp(ax/v)) is composed 
of solutions of the steady-state equation au - Vu = 0 and the method 
described above can be used to directly provide a difference scheme for 
the time-indpendent problem. The result is described by the two equations 
(II. 3a) and (II. 3c) in which the term is set equal to zero. This 

same basis can be used to construct a Green’s function on each overlapping 
subinterval having its singularity at x = x^, 

x^ ^ < x^ + ^, a technique which leads to a positive definite tridiagonal 

difference scheme for Sturm-Liouville problems as shown in [ 9 ]. An equiva- 
lent point of view has been independently developed and applied to similar 
singular perturbation problems which arise from steady-state problems (for 
example, c.f. [1], [3], [7], [10]). In this sense the methods described in 
this paper appear to provide the appropriate extention of such Green’s 
function techniques to time-dependent problems. 
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For the convective-diffusion equation a polynomial approximation 
basis also results by taking 

(f> x = <}>(0,a) = 1 

<f> 2 = 4>g(0,a) = x - at 

<f> 3 = 2 <j>gg(0,a) = i(x- at) 2 + V t. 

The difference scheme which is the consequence may be obtained from (II. 3) 
by setting p = q = 0. In view of earlier remarks this basis can be ex- 
pected to provide an accurate approximation only when the cell Reynolds 
number | 0 [ is small. 
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Conc luding Remarks 

This paper has described related implicit difference schemes for 
treating hyperbolic systems of equations and the scalar convective 
diffusion equation both of which share a common approximation rationale 
as well as a common solution technique. Numerical evidence indicates that 
both schemes can be employed to treat nonlinear problems. The accuracy 
of approximation to the dissipative hyperbolic problem is proportional 
to an artificial diffusion parameter a while the approximation to the 
convective-diffusion is second order accurate and is independent of the 
value of the local cell Reynolds number. For both schemes conventional 
energy estimates are available when the coefficients are constant. 

One important limitation of this study is its restriction to one 
dimensional problems; in another paper we will indicate how multi-dimensional 
problems may also be treated. 


Acknowledgment 

We wish to thank S. Wornom for valuable assistance with early numerical 
experiments and to V. Pereyra and J. Flaherty for helpful assistance about 
numerical methods underlying the solution of algebraic two-point boundary 
value problems. 



- 20 - 


References 

[1] A. E. Berger, J. M. Solomon, and M. Ciment, UnifioAmti/ accurate 
cU.66eA.ence. methods 6 0fl a singular perturbation pAoblem, Proc. Inti. 

Conference on Boundary and Interior Layers, Computational and Asymp- 
totic Methods, J. J. H. Miller, ed. , Boole Press, Dublin, 1980, pp. 

14-28. 

|23 C. deBoor and R. Weiss, Solveblok : A package 6 0fl solving almost 

block diagonal lineaA systems, with applications to spline approximation 
and numerical solution 06 ordinary di66crential equations, mrc Technical 
Report No. 16-75, University of Wisconsin Press, Madison, WI. 

[33 T. M. El-Mistikawy and M. J. Werle, Numerical method 6or boundary 
layer with flowing -the exponential box scheme, aim J., Vol, 16, 1978, 
pp, 749-751, 

[4] J. E. Flaherty and w. Mathon, Collocation with polynomial and 
tension splines 6 °r singularly perturbed boundary value problems, SIAM J. 
on Sci. and Stat, Computing, Vol. 1, No, 2, June I960, pp. 260-289. 

[53 H. B. Keller, Numerical methods in boundary layer theory. Annual 
Rev. Fluid Mech. , Vol. 10, 1978, pp. 417-433. 

[6] H. B. Keller, Numerical solution 06 two point boundary value problems, 

Regional Conference Series in Applied Mathematics, 24, SIAM, 1976. t 

[73 A. M. II' in, Vi66erencing scheme 6°r a di66erential equation with a 
small parameter a66ecting the highest derivative, Mat. Zametki, Vol. 6, 

1969, pp. 237-248, Math. Notes, Vol. 6, pp. 596-602. 



- 21 - 


[8] M. Lentinl and V. Pereyra, An adaptive, iinite diHenence SoZveA ion 
nontineoA two-point boundary pnobZems with mild boundary Zayens, SIAM J. 
Numer. Anal., Vol. 14, 1977, pp. 91-111. 

[9] M. e. Rose, finite. diiienence schemes Ion diiienentiaZ equations, 
Math. Comput., Vol. 18, No. 86, 1964, pp. 179-195. 

[10] M. e. Rose, Weak eZrnent appnoximation to eZZiptie diiienentiaZ 
equations, Numer. Math., Vol. 24, 1975, pp. 185-204. 

[11] M. E. Rose, A tneatment oi the wave equation and the Cauchy-Ziemann 
equations, to appear in SIAM J. Numer. Anal. 

[12] b. Wendroff, On centened diiienence equations ion hypenboZie systems, 
J. Soc. Indust. Appl. Math., Vol. 8, No. 3, 1960, pp. 549-555. 

[13] B. Wendroff, The stnuctune oi eeAtain iinite diiienence schemes, 

SIAM Review, Vol. 3, No. 3, 1961, pp. 237-242. 



Figure 1 


- 22 - 


aX=>/100 
\= 1 .04 


.40p 
.38 - 
.38 - 
• 34 — 


(la) 


•28 h 

.26 U 

-f 

•22 H- 
•2Ct 


■30 -3S .40 -4S .50 


*x =:/ :cc 

x- : .30 


• 40 p 
•38p 
■ 36*— 
.34 — 



.28 1- 
• 26i- 


(lb) 


.24 — 

•22 — 

•W U 
0 




*X=l/!00 

X= .70 


(lc) 


.40 r 

.38 — 
• 361— 
•34i- 
.32 — 



• 26 i— 


.26 j— 

.24 f- 

I 


•20 L 
0 



‘.2S 

X 


• 30 .3S -40 ’ -45 


Density profile of a shock calculation showing the influence 
of the CFL number upon pre- and post- oscillations (Scheme 1.6). 



Figure 2 
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. A comparison of the numerical and exact solution of a Riemann 
problem; the shock speed exhibits a relative error of 3% 
(Scheme 1.6). - 
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A comparison of boundary layer profiles for various cell Reynolds 
numbers . 
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Solutions of Burger’s equation 
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. Steady-state solution of Burger's equation. 
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